## get sra repository directory, then convert them to fastq file, then map it to reference genome to BAM file
## author: Yaping Liu  lyping1986@gmail.com
my $sra_dir = $ARGV[0];
my $prefix = $ARGV[1];
my $SRA_TOOLS = "/export/uec-gs1/laird/shared/research/yaping/software/sratoolkit.2.1.10-centos_linux64/bin/fastq-dump";
my $BSMAP = "/home/uec-00/shared/production/software/bsmap/default/bsmap";
my $SAMTOOLS = "/home/uec-00/shared/production/software/samtools/samtools";
my $PICARD = "/home/uec-00/shared/production/software/picard/default/";
my $ref = "/export/uec-gs1/laird/shared/research/yaping/data/genome_data/genome/hg19/hg19_rCRSchrm.fa";
my $mem="8g";
my $cores=6;
my $tmp_dir="/export/uec-gs1/laird/shared/research/yaping/data/tmp/";

my $header = "#PBS -q lairdprio\n";
$header .= "#PBS -l walltime=300:00:00,mem=$mem,nodes=1:ppn=$cores,arch=x86_64\n";
$header .= "#PBS -N sra_to_fastq_to_bam\n";
#$header .= "#PBS -W depend=afterany:1652168.hpc-pbs\n";
$header .= "#PBS -j oe\n";
$header .= "cd \$PBS_O_WORKDIR\n";	

opendir(DH,"$sra_dir") or die;
foreach my $file(readdir(DH)){

	if($file =~ /\S+.sra$/){

		$header = $header . "$SRA_TOOLS -B -F --split-3 $sra_dir"."$file -O $sra_dir \n";

	}
}
close(DH);

opendir(DH,"$sra_dir") or die;
foreach my $file(readdir(DH)){

	if($file =~ /\S+_1.fastq$/){

		$header = $header . "cat $sra_dir"."$file >> $sra_dir"."$prefix"."_1.fastq \n";
		
	}
	if($file =~ /\S+_2.fastq$/){

		$header = $header . "cat $sra_dir"."$file >> $sra_dir"."$prefix"."_2.fastq \n";

	}
}
close(DH);

my $read1 = $sra_dir.$prefix."_1.fastq";
my $read2 = $sra_dir.$prefix."_2.fastq";
my $outfileSAM = $sra_dir.$prefix."bsmap.hg19_rCRSchrm.sam";
my $outfileBAM = $outfileSAM;
$outfileBAM =~ s/\.sam/\.bam/;
$header .= "$BSMAP -a $read1 -d $ref -o $outfileSAM -z 33 -p $cores -b $read2 -s 16 -v 10 -q 2 \n";
$header .= "$SAMTOOLS view -bh -S $outfileSAM > $outfileBAM \n";
$header .= "$SAMTOOLS index $outfileBAM \n";
my $outfileBAMWithHeader = $outfileBAM;
$outfileBAMWithHeader =~ s/\.bam/\.header\.bam/;
$header .= "java -Xmx$mem -jar $PICARD"."AddOrReplaceReadGroups.jar I=$outfileBAM O=$outfileBAMWithHeader ";
$header .= "ID=$prefix LB=$prefix PL=illumina PU=run SM=$prefix CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT SORT_ORDER=coordinate ";
$header .= "TMP_DIR=$tmp_dir \n";

my $outfileName = "submitPBS";
open(OUT, ">$outfileName") or die "can not open file:$!";
print OUT $header;
#print $header;
close OUT;
system("qsub $outfileName");
print "$outfileName is submitted\n";
system("rm $outfileName");
	